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Abstract. Numerical simulations of turbulent stratified convection are used to study models with approximately the same 
convective flux, but different radiative fluxes. As the radiative flux is decreased, for constant convective flux: the entropy jump 
at the top of the convection zone becomes steeper, the temperature fluctuations increase and the velocity fluctuations decrease 
in magnitude, and the distance that low entropy fluid from the surface can penetrate increases. Velocity and temperature 
fluctuations follow mixing length scaling laws. 
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1. Introduction 

In modeling stellar convection it is important to make the 
models resemble the stars as much as possible. A major dif- 
ficulty in producing realistic simulations of deep stellar con- 
vection is the large ratio between thermal and dynamical time 
scales. This is because in dynamical units the solar thermal 
flux (Fq = 7 x 10 10 erg cm -2 s _1 ) is very small, 

; ^|~4xio- n , (i) 

pet 

at the bottom of the convection zone (where p = 0.2 g cm~ 3 
andc s = 200 km/s). At the surface this ratiois of order 10 _1 , 
but already two megameters below the surface the ratio is 
10~ 3 . This ratio is basically equal to the ratio of the dynami- 
cal to thermal time scales. 

For a common class of models studied by many work- 
ers (we will refer to them as polytropic models), most of the 
energy is carried by radiation instead of convection. This cor- 
responds to the polytropic index m of the associated hydro- 
thermal equilibrium solution having the value m = 1. In this 
case flow characteristics are different from those in the con- 
vection dominated transport regime, even for large Rayleigh 
numbers. In the range —1 < m <C 1, however, most of the 
flux is carried by convection. 

In order for simulations of convection to be representa- 
tive of late type stars the ratio of radiative to total flux must 
be small in the convection zone. This has prompted Chan & 



Sofia (1986, 1989) to study models where the radiative flux 
is removed entirely and is replaced by a subgrid scale flux. 
Unlike the radiative flux, which is proportional to the tem- 
perature gradient, the subgrid scale or small-scale eddy flux 
is proportional to the entropy gradient (e.g. Riidiger 1989). 
Others have chosen to model the surface layers and the gran- 
ulation (Nordlund 1982; Stein & Nordlund 1989, 1998; Kim 
& Chan 1998, Robinson et al. 2003, Vogler et al. 2005). In 
those models radiation is only important near the surface lay- 
ers and practically absent beneath the surface, although dif- 
fusive energy flux is still necessary for numerical simulations 
to be stable. 

The aim of this paper is to explore the effect of varying 
the diffusive radiative flux while keeping the convective flux 
in the convection zone approximately constant. From mixing 
length arguments one would expect that for negligible radia- 
tive flux the turbulent velocities, temperature fluctuations and 
other dynamical aspects of the convection only depend on the 
magnitude of the convective flux. However, this approxima- 
tion breaks down once the radiative flux is no longer very 
small compared to the convective flux. 

We consider models using piecewise polytropic layers. 
Such models have been widely studied and they possess the 
advantage that the properties of stable and unstable layers are 
easier to manipulate than in models with the more realistic 
Kramers' opacities, for example, where the radiative diffu- 
sivity depends on density and temperature rather than just on 



©2005 WILEY- VCH Verlag GmbH & Co. KGaA, Weinheim 



2.2 Radiative conductivity 



2 THE MODEL 



depth. However, we also include a comparison of realistic so- 
lar models with the same convective flux but varying diffusive 
subgrid scale energy flux in the interior and hence varying 
(turbulent) Prandtl number. 

There are two main results of our investigation. First, the 
properties of the convection depend on the ratio of radiative 
to convective flux when the radiative flux is not negligible. 
Increasing radiative energy diffusion reduces the tempera- 
ture fluctuations which requires larger velocities to carry the 
same convective flux. Increasing radiative energy diffusion 
also raises the entropy of the downdrafts and inhibits their 
descent. Second, as the resolution increases the dependence 
of convective properties on the Prandtl number decreases. 
However, at low resolutions the dependence is different when 
the viscosity is altered than when the radiative conductiv- 
ity is altered. Increasing the radiative conductivity to lower 
the Prandtl number has the effects described above. Decreas- 
ing the viscosity to lower the Prandtl number (increasing the 
Reynolds number) enhances small scale structures and in- 
creases the turbulence. 

We first discuss the setup of our model, the equations that 
are solved and how we change the fractional radiative flux in 
the convection zone while leaving everything else unchanged 
(Sec. |3- The polytropic model results are given in Sects [3] 
and |4] the realistic solar simulation results are presented in 
Sec. [5] and our conclusions are given in Sec.|6] 

2. The model 
2.1. Equations 

In our polytropic models we solve the continuity, momentum, 
and energy equations in non-conservative form, 
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strain tensor, v = const is the kinematic viscosity. We assume 
a perfect gas, so the pressure is given by 

V = (7 - l)pe, (5) 

where p is density, e = c v T is the internal energy, T is tem- 
perature, and c v is the specific heat at constant volume. 

Below the photosphere the radiative diffusion approxima- 
tion is valid so the vertical component of the radiative flux is 

F rad = KdT/dz, (6) 

where K is the radiative conductivity, T is temperature, 
and z is depth (increasing downward). We impose a cool- 
ing layer at the surface where r e (z) is a cooling time, and 
eo is the reference value of the internal energy e at the top 
of the layer. The value of eo is quantified by a parameter 
£0 = (7 — l)eo/(gd) = Hp° v ^ jd, which is the pressure 
scale height at the top of the box divided by the depth of 



the unstable layer, d. The pressure scale height is, H p = 
(dlnp/dz)^ 1 = TLT/(pg), where TZ is the universal gas 
constant and p the mean molecular weight. 

We adopt the basic setup of the model of Brandenburg et 
al. (1996, hereafter BJNRST; see also Hurlburt et al. 1986) 
where the vertical profile of JC(z) = K/c v is assumed such 
that K, is piecewise constant in three different layers, K = K\ 
in z\ < z < z%, K, — Ki in zi < z < Z3, and K, = JC3 in 
Z3 < z < Z4. In our case z\ — —0.15, Z2 — 0, 2:3 = 1, and 
z A = 2.8. 

2.2. Radiative conductivity 

In this paper we want to study the effects of varying the ra- 
diative flux in the convection zone, Z2 < z < 2:3, so we have 
to vary the value of /C2. In the following we discuss the rela- 
tion between IC2 and the anticipated radiative flux F la d in the 
convection zone, which will be a good approximation to the 
actual radiative flux obtained by solving Eqs. (|2ji-@. 

In astrophysics the radiative flux is often written in the 
form 

V = ■ z= 4- V. (7) 
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Here 



V = d\nT/d\np (8) 

characterizes the temperature stratification. We can define a 
radiative gradient, V ra d> that would be necessary to transport 
all the flux radiatively, 
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In the deep radiative interior of the sun we have F ra d = F to t, 
and therefore V = V ra d- In most of the convection zone the 
actual stratification is close to adiabatic, so V w V a d = 1 — 
7 _1 < V rad , since F rad < F tot . 

To set up a polytropic model it is customary to specify JC 
in terms of the polytropic index m (e.g., Hurlburt & Toomre 
1984), instead of V ra d- If all the flux were carried by radia- 
tion we would have p(z) ~ T m+1 , and so V ra d = (m+l) -1 . 
The value of the radiative conductivity is expressed in terms 
of the polytropic index from Eq. as 
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The three values of (i = 1, 2, 3) are given in terms of poly- 
tropic indices via the above equation. In all cases reported 
below we use TO3 = 3 and vary the value of 777,2. Near the top 
there is cooling (r e 7^ in z\ < z < Z2, and r e — > for 
z > zi). Therefore, almost all the flux in this layer is trans- 
ported by the corresponding cooling flux and the diffusive 
radiative flux in the uppermost layer is insignificant. Hence, 
the value of m\ does not affect the stratification. In most of 
the cases we put m\ — —0.9. 

The fractional radiative flux in the convection zone is 

1 -,(H) 



-Frad Vad / 1 1 . 
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where we have assumed V f» V a d- In all those models where 
777 = 1 (e.g. Hurlburt & Toomre 1984, 1986; Cattaneo et 
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2.4 The initial condition 



Table 1. Fractional radiative and convective fluxes for a few 
values of to, as obtained from Eq. (fTO . 



m 


-Frad/ftot 


-Fconv / Fiat 


1 


80% 


20% 





40% 


60% 


-0.5 


20% 


80% 


-0.8 


8% 


92% 


-0.9 


4% 


96% 



al. 1990, 1991; Brandenburg et al. 1990, 1996; Brummell et 
al. 1996, Ossendrijver et al. 2002, Ziegler & Riidiger 2003, 
Kapyla et al. 2004) and 7 = 5/3, the fractional radiative 
flux is as large as 80%. In order to study polytropic models 
with F ra( j <C F to t we must approach the limit m — > — 1 . In 
this paper we calculate a series of such models starting with 
to = 1 (the standard case considered in many papers) down 
to to — —0.9. As the value of to is lowered, smaller fractional 
radiative fluxes are obtained. In the upper layers of the solar 
convection zone this ratio is close to the ratio of the mean free 
path of the photons to the pressure scale height and can be as 
small as 10~ 5 . In deeper layers of the solar convection zone 
the radiative flux becomes progressively more important and 
reaches 50% of the total flux at 0.75 solar radii. 

It may seem unphysical to talk about negative values of 
to, because the density would increase with height, but such 
a 'top-heavy' arrangement only means that the stratification 
is then very unstable. Of course, the pressure still decreases 
with height, because m + 1 is positive. Negative values of to 
(but with m > —1) result from a very poor radiative conduc- 
tivity, which is indeed quite common in the outer layers of all 
late-type stars. 

It should be emphasized that to only characterizes the as- 
sociated thermal equilibrium hydrostatic solution, which is 
of course unstable for m < m a( \ = (7 — l) -1 = 3/2. In 
that case convection develops, making the stratification close 
to adiabatic. Thus, the effective value of m will then always 
be close to the marginal value m ad . The significance of the 
to used here (as opposed to the effective m) is that it gives a 
nondimensional measure of the fractional radiative flux. 

We carry out a parameter survey by varying the values of 
to 2 and hence /C e , which corresponds to varying V rac i in the 
convection zone and the total flux F to t . In practice we fix the 
value of F conv and calculate the total flux F from F conv for 
a given value of to-2. In Table^we give the fractional fluxes 
for some values of to, using Eq. ( II 11 1. Here we have assumed 
-Ftot = frad + ^conv, i-e. we have neglected kinetic and vis- 
cous fluxes which are of course included in the simulations. 



The flux is then expressed in terms of the non-dimensional 
quantity 

-ftot 
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which is a measure of the convective flux. (The actual con- 
vective flux, -Fc 0nv , is of course not a constant, but varies with 
height, and reaches a maximum of around JF conv .) 

The value of v is taken to be as small as possible for a 
given mesh resolution. For 50 3 we can typically use i/ = 5x 
10~ 3 (in units of yj gd 3 ). For a resolution of 200 3 we were 
able to lower the viscosity to v — 2.4 x 10~ 4 . 

The quantity T determines the 'mean' thermal diffusivity, 
X3 = W(7 A)), via 



where v|. ad is the radiative gradient in layer i, (Note that in 
this approach the actual thermal diffusivity decreases with 
depth in each of the three layers.) The nondimensional flux 
can also be related to the Rayleigh number, 

gd 4 fds\ 

Ra=^- — , (15) 

v X2 c p \dzj 

which characterizes the degree of instability of the hydro- 
static solution in the middle of the unstable layer, i = 2. Here, 

d fds\ 1-Vad/Vj 
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is the normalized entropy gradient of the associated hydro- 
static solution for the same JC(z) profile (see BJNRST). We 
can then express Ra in terms of T as 
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where Pr = v/\2 is the Prandtl number. In the astrophysi 

(2) 

cally interesting limit, VLj — > 00, we have 
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Ra -> Ra* = 2Pr 
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so large Rayleigh numbers correspond to small normalized 
fluxes and hence small velocities. This is at first glance some- 
what surprising, because large Rayleigh numbers are nor- 
mally associated with more vigorous convection and there- 
fore large velocities. However, while the velocity decreases 
in absolute units, it does increase relative to viscosity and ra- 
diative diffusivity, i.e. the Reynolds and Peclet numbers do 
indeed increase. 



2.4. The initial condition 



2.3. Nondimensional quantities 

Nondimensional quantities are adopted by measuring length 
in units of d, time in units of \fd]~g and p in units of its initial 
value, po, at z = Z3, i.e. at the bottom of the convection zone. 
In all cases we use a box size of L x = L y = 4 with Z4 = 2.8. 



The ratio of radiative to convective flux can easily be de- 
creased within the framework of standard polytropic models, 
provided the polytropic index m is close to —1. In that case, 
however, it is a bad idea to use the corresponding polytropic 
solution as the initial condition, because it is extremely far 
away from the final solution and would take a very long time 
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if 
if 
if 



Z\ < Z < Z 2 
Z-2 < Z < Z 3 
Z3 < Z < Z4 



(20) 
(21) 

(22) 



to relax to the final state. Instead, we use a simplified mixing 
length model with an empirically determined free parameter 
such that the entropy at the bottom of the convection zone and 
the radiative interior beneath is close to that finally obtained 
in the actual simulation. 1 

An initial condition that yields a solution on almost the 
right adiabat is obtained by assuming that the entropy is not 
constant, but that the negative entropy gradient is a function 
of the convective heat flux. According to mixing length the- 
ory the superadiabatic gradient in the convection zone is 

V-V ad = fc(F conv /pc s 3 ) 2 / 3 , (19) 

where A: is a nondimensional coefficient (related to the mix- 
ing length alpha parameter), which we determined empiri- 
cally for our models to be k = 1.5. 

In terms of e and In p, our initial condition can then be 
written as 

del dz = TO a( j3 V, 

dlnp/dz = m ad g (1 - V)/e, 

where 



V ad + k (F conv /pc 3 s ) 2 / 3 

and varies smoothly between the three layers. The initial 
stratification is obtained by integrating Eqs i20\ and (12 li 
from z — z\ downwards to Z4. The starting value of p at 
z = z\ is adjusted iteratively such that pizz) = po = 1, 
i.e. the density at the bottom of the convection zone is unity. 
This is especially useful if the extent of the box is changed, 
because this would otherwise affect the density in all layers. 

In Fig.^we compare the entropy obtained from the initial 
condition as derived above with that obtained from the ac- 
tual simulations. We also compare with the entropy derived 
from the associated polytropic hydrostatic thermal equilib- 
rium solution to show just how far away from the final state 
that solution would be. Finally, we also compare with a solu- 
tion where the entropy was assumed constant within the con- 
vection zone, which is better than the associated hydrostatic 
thermal equilibrium, but still with the wrong entropy at the 
bottom of the convection zone. 

As in previous work the initial velocity is random. This 
is realized by a superposition of randomly located spherical 
blobs of radius 0.1, where the velocity points in random di- 
rections. 

2.5. Boundary conditions 

In the horizontal directions we adopt periodic boundary con- 
ditions and at the top and bottom we adopt impenetrable, 
stress-free boundary conditions, i.e. 

du x du v 

— — = — — = u z = at z = zi,za. (23) 
oz oz 



1 The reason why we chose to specify our simulations in terms of 
polytropic indices is mainly that it makes good contact with previ- 
ous approaches using polytropic solutions. It is important to realize, 
however, that specifying the value of m in the convection zone is 
really just a way of specifying the nondimensional radiative con- 
ductivity. 
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Fig. 1. Entropy in the final and initial states, compared with 
the corresponding polytropic model, as well as a polytropic 
model with constant entropy within the convection zone 
(CZ). T com = 0.01. 

Table 2. Parameters for a model with J- conv = 0.01, £0 = 
0.2, i/ = 6x 10~ 3 and different values of m. 



m 


J tot 




Pr 


Ra 


1 


0.0500 


0.0400 


0.15 


9.3 x 10^ 





0.0167 


0.0067 


0.9 


2.1 x 10 4 


-0.5 


0.0125 


0.0025 


2.4 


8.9 x 10 4 


-0.8 


0.0109 


0.0009 


6.9 


3.3 x 10 5 


-0.9 


0.0104 


0.0004 


14 


7.4 x 10 5 



The boundaries are sufficiently far away so the boundary con- 
ditions do not significantly affect the flow properties in the 
convection zone proper, which is the layer z 2 < z < z 3 . 
Because of the cooling term in Eq. the value of e at the 
top is always close to the reference value erj which, in turn, 
is proportional to £0 and to the pressure scale height at the 
top. Thus, decreasing the value of £0 leads to stronger driving 
of convection and to stronger stratification in the top layer. 
However, fo cannot be decreased too much for a given nu- 
merical resolution, and so we chose £0 = 0.2 in all cases. In 
the top layer z\ < z < z 2 the inverse cooling time r~ x is 10 
and goes smoothly to zero for z > z 2 . 

3. Results 

3.1. Models 

We have carried out a series of calculations all with the same 
value of J-conv an d varying values of m 2 (hereafter referred 
to as m). In a first series of models we used c=6x 10~ 3 . 
The corresponding values of Prandtl and Rayleigh numbers 
are given in Table|2] 

Note that the Prandtl number is no longer small, except 
in the case m = 1. However, for m — 1 and J- COIW = 0.01 
the Rayleigh number is already so small that the instability 
to convection is suppressed. Therefore we have also studied 
another series of models where the convective flux is reduced 
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3.2 En tropy s tra tiUca tion 



Table 3. Parameters for a model with J- conv = 0.002 and 
v = 5 x 10~ 3 and different values of m. 

m Ttot -T^ad Pr Ra 

"LO 0.01000 0.00800 06 5.6 x 10 J 
0.0 0.00333 0.00133 3.8 1.3 x 10 5 



Table 4. Parameters for a model with J- collv = 0.01 and 
different values of m. v and \ 2 are given in units of 10~ 4 , i.e. 
i^_4 = v/l0~ A and \-4 = X2/10 -4 - The mesh resolution is 
also given. The asterisk in the last column indicates that the 
duration of the run was short and not yet fully relaxed, so the 
data may not be reliable. 



m 


V-4 


X-4 


Pr 


Ra 


resol. 




k T 


-0.5 


60 


25 


2.4 


9 x 10 4 


50 a 


0.34 


1.13 


-0.5 


12 


25 


0.5 


4 x 10 5 


100 3 


0.46 


1.07 


-0.9 


60 


4 


14 


7 x 10 5 


100 3 


0.29 


1.25 


-0.9 


12 


4 


2.9 


4 x 10 B 


100 3 


0.41 


1.14 


-0.9 


2.4 


4 


0.6 


2 x 10 7 


200 3 


0.48 


1.14 


-.98 


12 


0.8 


15 


2 x 10 7 


100 3 


0.38 


1.22 


-.98 


2.4 


0.8 


3.0 


1 x 10 8 


200 3 


0.43 


1.16 



by a factor of 5; see Table [3] Here however we have only 
calculated the cases m = and 1, because otherwise the ra- 
diative diffusivity would become so small that the resolution 
would be insufficient. 

Finally, we consider a series of models with varying 
Prandtl number and constant convective flux. Lowering the 
radiative flux automatically increases the Prandtl number. In 
the sun the Prandtl number is less than one, which was no 
longer the case in models with low radiative flux. We there- 
fore also studied the effects of lowering the Prandtl number 
by lowering the viscosity. In Table 0] we give the parameters 
for a series of models where Pr is varied by changing both v 
and X2- The question is whether or not certain properties of 
convection continue to depend on Prandtl number as the tur- 
bulence becomes more vigorous (small viscosity and radia- 
tive diffusivity). For instance, we expect that for sufficiently 
large Reynolds number the large scale flow properties will 
no longer depend on the microscopic values of viscosity and 
radiative diffusivity. 

3.2. Entropy stratification 

As the radiative flux is lowered, the mean entropy in the con- 
vection zone becomes more nearly constant and the supera- 
diabatic gradient at the top becomes steeper; see Fig.El The 
mean entropy beneath the convection zone is only slightly af- 
fected. 

As the radiative flux decreases the minimum entropy in 
the interior of the convection zone decreases (Fig. and 
Fig-EJ and approaches the minimum value that occurs at the 
top of the convection zone. These low entropy elements are 
only produced at the top where cooling is important. As the 
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Fig. 2. Comparison of the horizontally averaged entropy 
stratification for different values of m and T com = 0.01. 
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Fig. 3. Entropy histograms for different values of m and 
J collv = 0.01 at three different values of z. 



diffusive energy exchange decreases at least some fluid ele- 
ments make it all the way from the top to the bottom with 
very little mixing or heating. This is clearly a consequence of 
the reduced radiative diffusivity in the cases with low radia- 
tive flux. In the case m = 1 the deviations from the median of 
the entropy are rather small; see Fig.E] where .F conv = 0.002 
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m=-0.9 



m=0 




m= -0.8 




m= -0.5 





771= 1 




0.0 0.5 1.0 1.5 2.0 2.5 



Fig. 5. Entropy histograms for different values of m and 
•^conv = 0.002. Note that the range on the ordinate is reduced 
compared to the previous figure. 



A 




Fig. 6. Horizontally averaged density stratification for to 
Fig. 4. Entropy histograms for different values of to and n n , ^- _ n m 

_ — u.y ana j- conv — u.ui. 

J~ conv — U.UI. 



instead of T com = 0.01, so the entropy drop at the surface is 
obviously smaller. 

As the radiative flux is lowered, the drop of entropy and 
temperature at the surface becomes more sudden. This is be- 
cause with less radiative diffusivity (i.e. with less radiative 
energy transfer) the thermal boundary layer at the top be- 
comes thinner. At the same time the pressure must fall off 
smoothly, because the pressure is primarily determined by 
approximate hydrostatic balance. This can cause a rather pro- 
nounced density inversion near the top, as is seen in Fig. [6] 
which is occasionally also seen in the more realistic solar 
simulations. 



3.3. Convective and kinetic flux profiles 

The mean convective flux profiles (Fig.Q begin to converge 
to one and the same profile as the radiative flux is lowered. 
The differences between to = —0.8 and —0.9 are small, sug- 
gesting that with to = —0.8 or —0.9 the convective prop- 
erties of the simulations (convective velocities and tempera- 
ture fluctuations) begin to converge. However, the convective 
flux is not constant in the convection zone. This is because of 
some contribution of the kinetic energy flux, which is plotted 
separately in Fig. [8] As we lower the convective flux (from 
Jconv = 0.01 to J-conv = 0.002), the convective velocities 
decrease and therefore also the kinetic flux. 

The depth of the overshoot layer is clearly reduced as we 
decrease the convective flux (Figs0and|8j. This is in quali- 
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3.4 Velocity and temperature fluctuations 



S 1.0 
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> 0.5 - 
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0.0 0.5 1.0 1.5 2.0 2.5 



Fig. 7. Comparison of the convective flux for different values 
of m and T conv = 0.01 or 0.002 (marked by additional + 
symbols). 



Fig. 9. Comparison of the normalized vertical velocity fluc- 
tuations for different values of m and J- CO nv = 0.01 or 0.002 
(marked by additional + symbols). The normalization factor 

• -r— 2/3 
IS j conv • 



®5 




-0.25 
-0.30 



< 




0.0 0.5 1.0 1.5 2.0 2.5 



Fig. 8. Comparison of the kinetic energy flux for different 
values of m and T com — 0.01 or 0.002 (marked by additional 
+ symbols). 



tative agreement with theories linking the depth of the over- 
shoot layer to the magnitude of the convective velocities and 
hence to the convective flux (e.g. Hurlburt et al. 1994; Singh 
et al. 1998). 

3.4. Velocity and temperature fluctuations 

According to mixing length theory, the vertical velocity vari- 
ance is proportional to the relative temperature fluctuation, 
so {u 2 z ) ~ (AT/T) gt, where I is the mixing length, and 
gl ~ cf. The magnitude of velocity and temperature fluctu- 
ations (denoted below by primes) is determined by the con- 
vective flux via 



= {(pu z )'(c p Ty)~{ P ){ 



,2x1/2 



c p AT. 



(24) 



This estimate implies that, apart from factors of order unity 
(to be determined below), 



AT 



(ul) 



2/3 



(25) 



ct \ pet 

In Figsl9landll0lwe show the normalized velocity and temper- 
ature fluctuations for different runs: the profiles are basically 



Fig. 10. Comparison of the normalized temperature fluctua- 
tions for different values of m and .F conv — 0.01 or 0.002 
(marked by additional + symbols). The normalization factor 



-2/3 



is T, 



similar. Except for the cases m — 1 and m = 0, the mag- 
nitude of the velocity fluctuations decreases and the magni- 
tude of the temperature fluctuations increases as the radiative 



Eq. (1251 are indeed of order unity. It turns out that this is not 
only valid globally, but also locally; see Fig.^] Here we have 
defined the coefficients 



cz 



kr = 



(AT/T) 



(26) 



(T conv /(p C 3)}2/3' (Fconv/(pc 3 )) 2/3^ 

where the averages are taken over the convection zone proper, 
i.e. in Z2 < z < z%. 

The dependence of these scaling relations J25b and i26\ 
on Prandtl number is shown as a plot of k u and fcy on Pr in 
Fig- El The dependence of fcy on Pr is more or less unique, 
independent of whether the change in Pr is accomplished by 
changing v or Xi- By contrast, the dependence of k u on Pr 
is not unique. For constant values of v, k u is nearly constant, 
while it increases with decreasing values of \i- We explain 
this result as follows. 

As Xi increases (Pr decreases), temperature fluctuations 
are smoothed out, thus decreasing AT/T. This would de- 
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Fig. 11. Vertical profiles of the normalized mean squared ver- 
tical velocity fluctuations and temperature fluctuations, com- 
pared with the normalized convective flux raised to the power 
2/3. Note the good agreement between the three curves 
within the convection zone proper. 

crease the convective flux (which is proportional to the prod- 
uct of temperature and vertical velocity fluctuations), but 
since the convective flux is kept constant this can only be 
achieved by increasing the vertical velocity fluctuations and 
hence k u . 

The fact that k u and kx remain Pr-dependent even at rea- 
sonably large Reynolds numbers is somewhat surprising. Our 
largest Reynolds number, based on the rms velocity in the 
convection zone proper and the depth of the unstable layer, 
is around 1000. It is possible that one needs to go to much 
larger values before k u and kx become independent of Pr. 

3.5. Morphology 

Fig.^]shows temperature (or e) in horizontal planes at var- 
ious depths. At the top of the unstable layer the tempera- 
ture displays a familiar granular pattern with cool downdraft 
lanes. As m approaches —0.9 the pattern becomes generally 
sharper and more complex and of smaller scale. At the bottom 
of the convection zone (z = 1) there are isolated cool down- 
drafts. For to = —0.9 the granular surface pattern still pre- 
vails, but this is connected with weak density stratification. 
As the downdrafts enter the lower overshoot layer (z = 1.5) 
they become warmer. This is now due to the exterior entropy 
being lower than that of the downdrafts, which carry entropy 
from the upper convection zone. 

When the convective flux is decreased (JFconv = 0.002; 
see Fig. I14l >. the temperature pattern becomes sharper again, 
but the structure and the typical number of cells remains 
about the same. On the other hand, if v is lowered the tem- 
perature becomes much more intermittent and of significantly 
smaller scale; see Fig. [131 

4. Subgrid scale models 

Our next step is to compare our findings with the case where 
the radiative flux is replaced by a Smagorinsky subgrid scale 
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Fig. 12. Dependence of k u and kx on Pr. Note that v is 
constant along solid lines and x 2 and hence to are constant 
along dotted lines, v and \2 are gi ven m units of 10~ 4 , i.e. 
v-4 = v/lQ-* and = i//10~ 4 . 




Fig. 13. Horizontal slices of temperature at 4 different levels 
for four different values of to. J- CO nv = 0.01. 



(SGS) energy diffusion inside the convection zone (see Chan 
& Sofia 1986, 1989), corresponding to the limiting situation 
m = — 1. For this purpose, we use the original code of Chan 
and Sofia. Since the formulation (solving the conservative 
form of the Navier-Stokes equations) and the scheme of this 
code (second order only) are different from that used in the 
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Fig. 14. Horizontal slices of temperature at 4 different levels 
for two different values of to. !F conv = 0.002. 




Fig. 15. Horizontal slices of temperature at 4 different levels 
for m = —0.9 and three different values of v. T CO m = 0.01. 
A resolution of 200 3 meshpoints was used for Pr=0.6, and 
100 3 forPr=3and 14. 



previous section, we make a comparison of the results of the 
two codes using the previous to = —0.9, v = 0.006 case 
(with a 50 3 mesh). In Fig. the dashed and dot-dashed 
curves represent results from the Nordlund & Stein (1990) 
code (I) and the Chan & Sofia code (II) respectively. The 
agreement is good. 

In the SGS case, the radiative conduction and cooling out- 
side the convective region are fixed in the same way as dis- 
cussed in Sec. |2] Inside the convection zone, the numerical 
stability of the energy equation is maintained by a subgrid 
scale diffusive flux of the form = —\tpTVs, where \t is 
the SGS diffusivity. The numerical stability of the momen- 
tum equation is maintained by a SGS kinematic viscosity 

v t = 0.32 AxAz (2S 2 ) 1/2 . (27) 

where Ax and Az are the horizontal and vertical grid widths 
respectively. The ratio between v t and \t is fixed throughout 
the convection zone. 

In the original code of Chan & Sofia (1986) the effective 
Prandtl number, Pr t = Vt/Xu was chosen to be 1/3. How- 
ever, in order to facilitate comparison with the models dis- 
cussed previously, where we varied the thermal diffusivity, 
we now consider three models with Pr t = 1/3, 1, and 3. All 
three models have the same total flux F to t = 0.01 and use a 
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Fig. 16. Comparison of SGS models (solid and dotted lines, 
m = — 1) with direct calculations (to = —0.9). The dashed 
and dash-dotted lines refer to Code I used in Sec. [5] (Nord- 
lund & Stein 1990) and Code II used in the present section 
(Chan & Sofia 1986). While the to = -1 and to = -0.9 
models are very similar in the mean density stratification and 
the scaled temperature fluctuations, there are significant dif- 
ferences in the mean entropy profile and the scaled vertical 
velocity fluctuations. 



50 3 mesh. In the SGS case the flow still shows the usual gran- 
ular pattern at the surface; see Fig.^] Smaller diffusivity xt 
produces thinner intergranular lanes and decreases the size 
of the smallest granular structures; it is consistent with the 
trend shown in the previous section. Compared to the con- 
stant v case with the same mesh size (Fig.ll3l>. the SGS pat- 
terns show more vortical features and indicate more vigorous 
turbulence. In the case with the largest Pr, however, the tem- 
perature field is getting a little noisy; the thermal diffusivity 
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1.2 



Fig. 17. Horizontal slice of temperature for m = —0.9 and 
SGS model using three different values of the Prandtl number 
Pr = vt/xt- -^conv = 0.01. All calculations here use a 50 3 
mesh. 



is already too small to give reliable results. In Fig.[H)]we plot 
the mean entropy, density, as well as velocity and temperature 
fluctuations for Pr t = 1/3 and 1. 

Since v t depends on the local strength of the turbulence, 
it is not uniform anymore. Its horizontal mean varies with 
depth and reaches a maximum near the top of the convection 
zone. The range of the mean values (0.0015 — 0.0030), how- 
ever, is quite limited and they are less than half of that used 
in the previous comparable case (m — —0.9, v — 0.006 and 
same mesh size). Relative to this case, the effective Reynolds 
numbers of the SGS models are thus considerably larger. 
The SGS diffusivity xt, on the other hand, (oc Pr^ 1 ap- 
proximately) varies by almost an order of magnitude across 
the different models and are larger than the x of the low- 
resolution to = —0.9 case (ratio w 1.1 — 10). The temperature 
fluctuations of the SGS models are thus smaller (difference 
f» 20 — 30%). To deliver a comparable amount of convective 
energy flux, the value of the vertical rms velocity for the SGS 
models is considerably higher than that of the m = —0.9 
model. This is accomplished by a somewhat steeper entropy 
gradient in the convection zone. The density stratification is 
similar. Given that the turbulence has become more vigorous, 
the kinetic energy flux is now increased by almost a factor 
of 2 to about 35% of the total flux. Part of that (~ 9%F tot ) 
is compensated by the subgrid scale convective flux, and the 
rest is balanced by the enhanced enthalpy flux. 

5. Solar simulations 

We have also made two simulations of near surface solar con- 
vection with Prandtl numbers Pr t = Vt/xt = 0.33 and 3.3, 
using realistic physics, with a low resolution of 63 3 (for de- 
tails see Nordlund & Stein 1990; Stein & Nordlund 1998). 
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Fig. 18. Total, enthalpy, radiative and kinetic energy fluxes 
for the solar simulations. Crosses are case Pr t = 0.3 and lines 
are case Pr t =3. When the diffusive energy flux is larger, the 
kinetic energy flux increases and the enthalpy flux decreases 
slightly. 
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Fig. 19. Average entropy vs. depth. (Here the Prandtl number 
is written as a.) Larger energy diffusion produces a slightly 
smaller entropy gradient just below the surface. 



Here the radiative flux vanishes below the surface but there is 
numerical energy diffusion with the above ratios to the mo- 
mentum diffusion. The vertical energy diffusion is propor- 
tional to the gradient of the energy minus the horizontally 
averaged mean energy, so it is similar to the entropy diffu- 
sion used in the SGS case. The diffusion coefficients are not 
constant, but have terms proportional to the sound speed, the 
magnitude of the velocity, and the compression. They are en- 
hanced where there are small scale velocity fluctuations and 
quenched in laminar regions by the ratio of the magnitudes 
of the third to the first derivatives of the velocity. These sim- 
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Fig. 22. Horizontal slices showing the temperature at depths 
Fig. 20. Average density vs. depth. Larger energy diffusion of 0.1, 0.5 and 2.0 Mm for two values of the Prandtl num- 



leads to a less extended atmosphere. 
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Fig. 21. Histogram of the entropy distribution at a depth of 
1.5 Mm. For larger diffusive energy transfer, smaller Prandtl 
number, the lowest entropy fluid is destroyed and the expo- 
nential distribution becomes steeper. 



ulations only extend to 2.5 Mm below the surface and the 
convective flux is controlled by specifying the entropy of the 
inflowing fluid at the bottom. The simulations were run for 2 
solar hours. 

Larger energy diffusion (i.e. smaller Prandtl number) pro- 
duces a slightly larger kinetic energy flux and a slightly 
smaller enthalpy flux, resulting in a 5% reduction in the net 
flux (Fig. II 8t . The energy transport switchover between ra- 
diative and convective occurs slightly deeper in the large en- 
ergy diffusion case. 

There are some other small alterations in the mean struc- 
ture: larger energy diffusion produces a less steep mean en- 



ber. Dark is low temperature and light is high temperature. 
Each panel is scaled independently. Smaller energy diffu- 
sion, larger Prandtl number, allows smaller scale temperature 
structures to exist. 
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Fig. 23. Horizontal slices showing the vertical velocity fluc- 
tuations at depths of 0.1, 0.5 and 2.0 Mm for two values of the 
Prandtl number. Dark is upward velocity and light is down- 
ward velocity. There is little difference in the scale of the ve- 
locity structures between these high and low energy diffusion 
cases. 



tropy gradient at the surface (Fig. \\9\ and a less extended 
atmosphere ("Fig.l20>. 

The low entropy fluid (which gives rise to the buoyancy 
work that drives the convection) is fluid that reaches the sur- 
face and radiates away its energy and entropy. Larger en- 
ergy diffusion destroys the lowest entropy fluid as it descends 
back into the interior, by heating it up. This leads to slightly 
steeper exponential decline in the entropy probability distri- 
bution function and a less extended low entropy tail to the 
distribution (Fig. l2U . 

Increasing the energy diffusion has a clear direct influ- 
ence on the temperature structures in these simulations, just 
as found when varying the diffusive radiative flux - larger 
energy diffusion produces larger more diffuse temperature 
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Vorticity (s ') 

Fig. 24. Histogram of the vorticity distribution at a depth of 
1.5 Mm. Varying the energy diffusion by a factor of ten has 
only a very slight effect on the velocity and vorticity. 

structures, smaller energy diffusion allows smaller, sharper 
temperature structures (Fig.l22>. 

The velocity and turbulence, on the other hand, are little 
affected by this factor of ten variation in the energy diffusion 
(Fig. l24l >. Changing the resolution and hence the viscous mo- 
mentum diffusion, however, has a profound affect on the tur- 
bulence (vorticity) and the velocity (Stein & Nordlund 1998). 
Less viscosity leads to greater turbulence, larger vorticity, as 
well as a velocity distribution extending to larger magnitudes 
in all directions but only in a small fraction of the volume. 

6. Summary 

There are a number of clear changes in convective proper- 
ties as the diffusive radiative flux is decreased while keeping 
the convective flux constant in a convection simulation. First, 
the entropy jump near the top becomes larger and steeper and 
the low entropy fluid produced by cooling at the surface pen- 
etrates farther through the convection zone leading to a fi- 
nite probability to find small regions with very low entropy 
near the bottom of the unstable layer. Second, the temper- 
ature fluctuations increase and the velocity fluctuations de- 
crease in such a way that their product, which is proportional 
to the convective flux, remains approximately constant. Also, 
the kinetic energy flux decreases. Third, the dynamics in the 
overshoot layer becomes somewhat more intermittent due to 
a few strong downdraft plumes. Finally, in all cases the veloc- 
ity and temperature fluctuations follow mixing length scaling 
laws; see Fig. II II 

The radiative flux really serves two different purposes: 
it transports heat vertically, and it keeps the model numeri- 
cally stable by diffusing energy fluctuations both horizontally 



and vertically. Since those two properties appear to be rea- 
sonably well decoupled from each other, one might separate 
them by having a small vertical radiative flux plus a subgrid 
scale diffusive flux that keeps the model stable. It is in prac- 
tice difficult to decouple the need for energy diffusion from 
the vertical diffusive heat transport, which one would like to 
keep small if one diffuses on the temperature. However, as 
discussed in the introduction, such a separation is possible if 
the diffusive flux is based on entropy or temperature fluctua- 
tions. Convective simulations with small radiative fluxes, as 
is appropriate for cool stars, would then be feasible. This was 
recently demonstrated by Miesch et al. (2000) and Brun et al. 
(2004) using simulations of fully spherical shells. 
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